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АВ5ТКАСТ 


Many meteorological forecast applications require the use of grids 
that have a high resolution in a particular area of interest, while 
allowing coarser resolution elsewhere. Conventional finite difference 
models often use nested grids to this end. In recent years, finite 
element models have been offered as an alternative. In this study, the 
two-dimensional advection equation with diffusion is defined over a 
rectangular domain. The Galerkin technique is applied to linear basis 
functions on triangular elements. The model is tested to determine the 
sensitivity of the forecast to various nodal geometries. Both equi- 
lateral and right triangular elements are tested. It is found that the 
equilateral arrangement consistently yields a superior forecast. Other 
tests are conducted in which the resolution is varied smoothly versus 
abruptly over the domain. The smoothly varying case gives results that 
are dramatically improved over the abruptly varying case. Among the 
conclusions is the fact that, for a given maximum resolution, the more 
Slowly and smoothly the element size is changed, the better the forecast 


obtained. 
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Т. INTRODUCTION 


Sub-Synoptic scale meteorological elements are usually not repre- 
sented or forecast well by coarse mesh numerical models. This situation 
may be improved by increasing the resolution of the grid. However, a 
uniform reduction of the grid mesh requires a significant increase in 
computational time and computer storage. Hence this is a limitation on 
the practical resolution, and therefore accuracy, for a uniform grid. 
The conventional solution is to nest grids so that the geographical area 
or meteorological feature of interest 1s covered by a fine mesh, and 
less interesting areas or features are left with a coarse mesh. The 
advantage of this technique is to increase accuracy in the desired area 
while keeping the computational storage and time requirements within 
reasonable limits. There are, however, two disadvantages to nested 
grids. First, the abrupt change in grid size at the boundaries of the 
fine mesh generates noise in the solution. Second, boundary conditions 
for the fine mesh must be interpolated from the coarse mesh. 

A better technique would allow the grid to vary smoothly and con- 
tinuously over the domain rather than in abrupt jumps. However, such a 
highly variable grid greatly complicates a finite difference numerical 
model. As an alternative to such a finite difference technique, consider 
the Finite Element Method (FEM). This method has long been used in 


mechanical engineering but has been adapted to meteorology only during 


Il 


the last decade. Pioneering work in meteorological applications was 
Mmm” Cullen (1973). Staniforth and Mitchell (1978) and Staniforth 
and Daley (1977) have devised more recent forecast models. The most 
recent FEM meteorological model at the Naval Postgraduate School was 
written by Kelly and Williams (1976) and is the precursor to this study. 

The great attractions of the FEM are the more accurate phase speeds 
[Haltiner and Williams (1980)], and its suitability to non-uniform grids. 
The method involves the division of the domain into a number of elements 
and then defining an equation for each node of each element. Thus, 
solution requires solving a large number of linear equations simul- 
taneously. However, there is no limitation on the size of the elements 
or the shape of the domain. The method is easily adapted to concentra- 
tion of small elements in areas of interest and assignment of larger 
elements to other areas. 

The purpose of this investigation is to determine the sensitivity of 
a FEM model to different grid geometries. The model used is a two- 
dimensional advection model with diffusion. It was deliberately chosen 
because of its simplicity. By keeping the model simple, all variations 
in the results can be attributed to differences in geometric and initial 
conditions rather than inaccuracies in the equation formulation. In 
addition to sensitivity studies, the model can be adapted as an opera- 
tional predictor of any scalar field such as cloudiness, vorticity, 


temperature or precipitation. 
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The Galerkin technique is employed to transform the differential 
equation for the dependent variable into a set of algebraic equations. 
Triangular elements are used with the area coordinate system. Later 
chapters will show that these techniques greatly simplify the 


computations. 
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The Galerkin technique will be explained by first an example, after 
Crandall (1956), and then the application to the FEM, after Kelly and 


Williams (1976). 


A. EXAMPLE 
Consider a simple differential equation, that which governs the cool- 


ing of an object 


шиг al Е = 0 
E. m ЕТО 
dt 


with the solution defined on the interval 0 < t < 1. 
The critical step is to select a trial family of approximate solutions. 
(The members of a trial family are often called basis functions.) For 


Simplicity, consider the second order trial~solutions 


^^ 2 
= t E F I- 
x T с, cot II-2 


^ 
where x is the approximate solution. 


The object is now to determine the coefficients c, and c,. Next, form 


1 2 
the residual, R(t) from II-l 
dx 
К = — + - 
3t х II-3 
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substituting II-2 into II-3 


Ber c, dmt) + Es ГЕТЕ t?) II-4 


Now we want to adjust сі апа с. so that II-4 stays close to zero on the 
interval O «€ t < 1. 
The above discussion applies to any weighted residual technique; 
however, Galerkin's technique has two requirements: 
1. Тһе wetghted averages of the residual must vanish over the interval. 
2. The weighting functions must be the same functions of t as were used 
to construct the trial family (in II-2). In this case those 
weighting functions are t and m 
That is, the weighted averages are formed by multiplying the residual 


by the weighting functions, integrating the product over the interval 


and setting the result equal to zero. 


F2 


! 
© 


ERO dt 


and II-5 


Ен 


|! 
O 


tem at 


ES 





Sepstituting II-4 into II-5 


1 i 
2 3 
t Rdt = | [t * c, (t*t ) + с, (26 +t )jdt 
о о 
1 5 D 
= = + = + — = - 
5 с с. 12 с. O II-6a 
and 
B ШІ 
ЕТЕ ele = m c(t +t) + c, t^ t^) Jat 
O O 
1 7 7 
=—+ > = = 
3 12 с, 18 с, 0 II-6b 


Solving II-6a and b for с, апа с, we find: 


с, ш -.9143, с. = 0.2857 


So that the approximate solution is 


leas E + 0.2857 t^ 


The exact solution to Equation II-1 is 


l6 


Figure 1 shows the plots of the exact solution and the approximate solu- 
tion. Note that this technique yields the best fit only on the desired 
interval O < t < 1 and that the solutions may diverge from the interval. 
In general, Galerkin's technique will reduce a partial differential 
equation to a family of N algebraic equations, where N is the number of 
basis functions. In this example Galerkin's technique has reduced II-1 


: ; : 2 
to two equations, II-6a and II-6b, in the two basis functions t and t. 


B. APPLICATION TO THE FEM 

The Finite Element Method divides the domain into discrete segments 
called elements with points (called nodes) arranged about the perimeter of 
the elements. The basis functions are then defined locally. These basis 
functions are usually low order polynomials that must be piece-wise 
continuous. A one-dimensional example is Figure 2 wherein the domain 
(x-axis) is divided into four elements (line segments) W through Z. All 
of the basis functions (A-C) are linear rising to a value of unity over 
each node (points 2 through 4) and are zero elsewhere., Now, if a variable 


S is defined over node 3, for example, then it may be approximated by 
= И + -7 
o (с = 


where Yor Y3 and А are the values of S at nodes 2, 3 and 4 respectively 


and V 


AU V 


T Ve are the basis functions A, B, and C respectively. 
Equation II-7 may be substituted into any equation in which S appears. 


An equation using II-7 can be written for each node, by multiplying the 


basic equation by the basis function and integrating over the domain. 
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Figure 1. Solutions to cooling object example. Approximate solution 
is the curve connecting the 'x' signs. Exact solution is 
the curve connecting the '+' signs. 
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Figure 2. One-dimensional example of linear basis functions. 


5 6 Y 
/ | 
х 
2 3 


Figure 3. A domain with eight nodes and eight triangular elements. 
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The end result is a system of N equations in N unknowns, where N is the 
number of nodes. Fortunately, in matrix form, these equations are often 
tridiagonal and may be solved with great economy. 

In two dimensions the basis functions are somewhat different, but the 
mathematics is identical. In Figure 3, the domain has been divided into 
eight triangular elements with eight nodes. Figure 4 shows the basis 
function (outlined in heavy black) at node 7 that must span all the 
elements about node 7. Note that the function has value of unity at 
node 7 and decreases linearly to zero at all surrounding nodes. Figure 5 
shows another basis function, but this one at node 5. Each node has a 
Similar basis function about it. 


The value of a variable S may be approximated at any node i by 


S = "S II-8 
E. 


J 


where j ranges over all the nodes connected to node i including i itself, 
Y; is the value of S at node j and Ws is the value, at the j-th node, 
of the basis function of the i-th node. 

As in the one-dimensional case, an equation like II-8 is written 
for each node, i, then multiplied by the basis function and integrated 
over the domain. The resultant equation set, in matrix form, contains 
an n by n square matrix where n is the number of nodes. Although this 
matrix is usually not tridiagonal, careful selection of a nodal numbering 
System would yield a matrix that is strongly banded and fairly easily 
solvable. Alternatively, with these sparse matricies, iterative methods 


converge rapidly. 
20 








Figure 4. The linear basis function associated with node seven. 





Figure 5. Another linear basis function on the same domain. 
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III. EQUATION FORMULATION 


The equation of interest to the investigation is the two-dimensional 


advection equation with second order diffusion: 


95 95 95 2 ОЖК 
SE + un Bu ss - К "8 = 0 111-1 


Define an inner product of two functions f(x,y) and g(x,y) as 


| [е == Erg. 
y x 


Following the Galerkin technique, form the residue and find the weighted 
averages, uSing the basis functions as weighting functions. Neglect the 


diffusion term at this stage. 


82 + u 92 + у ЫН V dx dy = 0 
ух 
or 
д 
«5, у> + ы 35, у tev Sw = 0 III-2 


where V is the weighting function. Although the domain of integration 
is over the entire range of x and y, it will soon be shown that V is zero 
everywhere except locally around a specific node. So that for any node, 


III-2 is non-zero only when within one grid increment of that node. 
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The variables will have the form 


5 = ү, У, 
5 2 

а= о Y, III-3 
J J] 
DEUS WV 
ge] 


where the repeated subscripts indicate summation over the range of the 
subscript. The coefficient is a function of time only and the basis 


Buucblon js a function of space only. That is 
б= ү, У. = } (EI MEN) 
Y} Yj j ,У 


etc. 
The critical step is the requirement that the weighting function in III-2 
be the same as the basis function in III-3. Substitution of III-3 into 
III-2 yields 
. 9У, 
ни? s 


ӘУ, 


J 
< 2 a i > — = 
se Bk Vk Y, 3 , л 0 ІІІ-4 


where the dot indicates a time derivative. 
Applying first the Galerkin technique and then the Gauss Divergence 


Theorem to the diffusion term yields 
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Куе ахау К ЇЇ * (VS)V, dxdy 
y x y X 


К, [= - VS + W,] ахау Tre: 
y x 


Kt vs. n dr - VV. * VS dxdy] 


where n is a unit vector normal to the domain and dr is the differential 
distance along the path of integration on the perimeter of the domain. 
95 
The contour integral in III-5 vanishes because = cancels out due to 
En. 95 
elec continuity o O because there is no flux across the channel 
walls at the north and south edges of the domain. Using III-3, the re- 


maining term on the r.h.s. of III-5 can be written 





ду, Ir ду, ду, 
-K, «WW, * VS? = = A 2 ш + MS : Y, mm III-6 


Combine III-4 and III-6 into the transformed advection equation as 


follows 
е or 22 
AV, > + Y 10 < — ‚> < E > 
ШЕ m i Yad kk oe” u > 5, Er 220 E 
9v. NE OV x 
t ‚> t <S PERS 22 = 
585003, y 
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Now using a centered difference in time gives the forecast equation 


ntl n-1 mS n n < J > п «y re. ү > 
(ү; к. EU 2At Y, ta, UN SSMO 8. ee 
A ov, ду, QUE 
сан m eee иг, Ї11-7 
ш к! Өх ' Ox ду ' ду 1j 


when n is the time step. In this model a” and Be are not forecast 

quantities, but are either specified for each time step or are constant. 

They represent the prescribed wind field which is advecting the scalar S. 
This equation is valid for each node i and as such is a matrix 


equation of the form 


[A] {x} = {b} 


where 


[A] = < 9 : Vi > with dimensions n by n 
Ф ar T . “ . 
{х} = ( E 2 - Ya i ) with dimensions 521 


{b} = The right-hand side of III-7 
with dimensions 125201 
апа п 15 the number of nodes. 
Note that all of the inner products are independent of time so that 
they need only be calculated once and placed in mass storage to be 


accesSed at each time step. 
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IV. EVALUATION OF INNER PRODUCTS 


Even though the choice of elements in this investigation is triangles 
with nodes at the vertices, this is by no means the only choice possible. 
The elements could be triangles with six nodes each; rectangles with six, 
eight or twelve nodes; or any other polygon or curved-sided element. More 
complex elements would allow higher orders of continuity.  Bathe and 
Wilson (1976) describe the technique of isoparametric interpolation as 
applied to the FEM. However, these more complex elements require even 
more complex formulations. The mathematics can get out of hand very 
quickly requiring many more terms and complicated integrations. Also 
Galerkin's method is not the only weighted residual technique available. 
Crandall (1956) describes the collocation, subdomain, and least squares 
methods as well. Lower order continuity is deemed to be acceptable in 
this investigation because the grid resolution is fine enough so that 
the change in variable is assumed to be linear between grid points. 

This assumption allows the selection of triangular elements and the 
Galerkin method, each of which will significantly simplify the mathematics. 

Figure 6 shows a typical element with a point p described in area 

coordinates. Lines are drawn from p to each of the vertices 2 Zr: 


j=1,2,3) dividing the element into the areas 2r (1:1:22:31: 


Define: E. Бо 


where A total element area 
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P (x,y) 
or 
P(L,,L,,L;) 


X1 ‚У; 


X2,Y2 


Figure 6. A typical element illustrating the area coordinate 
system. 





Figure 7. An element with part of a basis function. 
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Now 


+ + = 
Li L, L 1 
= = т 
х LX. Lx. Lx, 
= + 
ЖЕГЕ“ 5272 5 54/5 
or, in matrix form 
ү 1 1 1 L 
l 
x = Xx х, х, L, 
Чи Уі Y > Y 3 E, 
Srpesolving for the L vector 
Бү 2A ру а, l 
L m 2A b a ІУ-1 
2 ~ 2A 2 2 2 
L 
3 2A р. an Y 
where a, = X4 - X, b} ECT E. 
aces X, = xX b, = y, = y 
2 L 3 (7112 3 1 
юэ. 5 537 3) 7 У 
So that any one of the equations in IV-1l is really 
j E 
EE yc — xr — = 
j ne on IV-2 


28 





But by the chain rule 


DN» 732 
х 9х ðL, 2A OL 

and IV-3 
eros e 


“Милл” ры алуы eee D 


where the repeated subscript indicates a summation from one to three. 
Consider Figure 7. The heavy black lines outline an element 
divided into area coordinates defining point P, and only A, is labeled. 
The hatched triangle (labeled у.) 1s that part of the basis function 
associated with node 2 that lies over the element. There are two more 
basis functions called Vi and V4 that are associated with nodes 1 and 
3 respectively. Both of these functions also have linear sections over 
the element. However, they have been omitted from Figure 7 to avoid 
clutter. The altitude h, is defined as unity. The altitude h is the 
value of V, at P. Recall L, = A/A. As P moves from node 2 to the 


opposite side of the element, h varies linearly from 1 to 0. Following 


the same path L. also varies linearly from 1 to 0. So that at any 


2 


point P on the element L, = h= V. In general, te = ЫЕ This is the 


first great advantage of triangular elements. 
So now, using IV-3 


V 
д і b у, b, 9У, 5, oV 


A 
Ox 2 oL, 2A oL. 2A oL. 


5, IL. b, ðL, b, 91. 


2А oL, 2A 91. 2А dL. 
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but 











IL. 
3L. = 1 ТЕ 1 = jJ 
J 
= 0 if i ž# j 
So 
9v. b. ду, а. 
5; ES and similarly am EE IV-4 


Note that these derivatives are dependent upon the individual elements 
and are constant with time. Therefore they need only be calculated once 
and be stored as part of the inner products. 

The second great advantage of triangular elements and linear basis 


functions is the general area integration formula 


m n min! 
Da L, aA = ас 50) 2А IV-5 


Or, rewritten as an example of an inner product: 





2 211! A 
« > = = 
= í > 51 = 30 


Now the inner products for equation III-7 may be rewritten with the 


aid of IV-4 and IV-5: 
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V. INITIAL AND BOUNDARY CONDITIONS 


Simple initial conditions are used for the scalar field S and the 
wind components u and v. The S field consists of a half sine wave in 
the y direction and multiple waves in the x direction. These multiple 
waves are the mechanism used to generate smaller scale features. The 
= term in the u equation causes the area of fastest flow to be coinci- 


dent with the area of highest spatial resolution. 





S=As E ) >) V-1 
u = U + B cos (ux ш) V-2 
L 2 


where 


A = arbitrary amplitude = 400 meters 


n = wave numbers SL oo 
L = channel length = 2400 km 
W = channel width = 2400 km 
U = mean zonal wind = 20 m/sec 


B= U x R; R is the ratio of the perturbation 
wind to the mean zonal wind 


Figures 8 and 9 depict the initial fields of S and u resvectively. 
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Figure 8. Initial field for the scalar S. 
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Figure 9. Initial u field with wind ratio R = 0.4. 
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Boundary conditions are also handled simply. Cyclic continuity is 
assured in the x-direction by selecting the proper node convectivity. 
Two techniques are used to control flow across the channel walls at 
the north and south sides. First, the exponent on the y-component of 
the initial wave is high enough to keep the value of S close to zero 
in the vicinity of the walls. Second, the solution method is iterative 
so that the values of the nodes one row in from the channel walls are 
Set equal to the nodes on the walls at each iteration. The result is 


that diffusion from the edges into the interior is disallowed. 
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VI. COMPUTATIONAL TECHNIQUES 


ee ЫЕ = ө 


A. NUMBERING SCHEMES 

As described in Chapter IV, the inner products may be easily 
evaluated on the individual element level. To facilitate this evalua- 
tion, a local numbering system is required for each element. An array 
called IEL, dimensioned N by 3, contains the identification number of 
the nodes on the three verticies for each of N elements. For each 
element these node numbers must be stored sequentially in a positive 
Sense, that is counterclockwise. The ncde with which to start number- 
ing, however, is arbitrary. 

In addition to the local numbering system, a "mass" matrix is 
needed. The mass matrix is a matrix of coefficients whose rows are 
the equations of the system to be solved. There is one equation (row) 
for each node. Each equation will have a term (column) for each node. 
A non-zero entry in the mass matrix indicates connectivity. That is, 
if entry (I,J) is non-zero then node I is connected to node J. Two 
nodes are connected if they are both vertices of the same triangular 
element. In this model, each node is connected to a minimum of two and 
a maximum of six or eight other nodes depending upon model version. 
The model uses an array called MAME dimensioned N by 7 (or 9) that con- 
tains the numbers of all nodes connected to any specific node, includ- 
ing itself, for each of N nodes. MAME will contain from 3 to 7 (or 9) 


non-zero entries in each row. 
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To save storage space, MAME is next compacted into a one dimensional 
vector NAME which contains only the non-zero entires of MAME. The 
utility vectors ISTART(N) and NUM(N) are also constructed in order to 
decode NAME. ISTART(N) contains the index in NAME of the starting 
position of the continuity of any node М. NUM(N) contains the length 
of the connectivity of node N within NAME. As an example of how to re- 
trieve the connectivity of a node, say node 14, from NAME: Let J - 
ISTART(14) and let K = NUM(14). Then NAME(J) through NAME(J+K-1) contain 
the connectivity for node 14. 

Consider the matrix equation [A] {x} = {b} from Chapter III. It 
should be clear that [A] is the mass matrix mentioned above. It is a 
very sparse matrix containing the inner products BL Mam for row j and 
column i. In this model [A] is compacted into the vector {H} using NAME 
as a lookup table. The result is that there is no longer a real matrix 


equation, but rather just the product of vectors 


ex = i} магт 


B. SOLUTION TECHNIQUE 

An advantage of economy of space has been gained, but VI-1 now re- 
quires an iterative solution. The solution chosen is a standard Gauss- 
Seidel technique which requires a reversal of direction of iteration 
every pass. An average of 16 passes are needed for each time step. The 


solution is considered to have converged to its final value when: 
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m every node i. 


with: x iteration number 


solution for node i in x 


E 
18 


More efficient solution schemes are available, but this one was 
selected because of its simplicity. Minimizing computation times is 
not of high priority to this investigation for two reasons: first, the 
recently installed IBM 3033 is a fast machine, and Secondly, since an 
advection model doesn't contain gravity waves, stability can be main- 
tained with very large time steps. As an example, a twelve hour 
forecast with a grid size of 200 km and a time step of one hour requires 
less than eight seconds of CPU time. Even the largest version of the 
model with longer forecast length, higher resolution and linkage to 


the plotting software takes less than eight minutes. 


ERZZROBERT FILTER 

Haltiner and Williams (1980) discuss the advantages of time averag- 
ing. Averaging operators act as low-pass filters that tend to remove 
high frequency waves while having little effect on long-period waves. 
Even very weak filters will remove high frequency noise if used at every 


time step. In this model, the choice of filters is that of Robert (1966). 


37 





First, assume that an average field ER 1! already exists at time 
level n-l as well as the unaveraged field 5, at time level n. Then 
the model uses its predictor equation to determine the solution vector 


{x}. For each node, an unaveraged predicted value is computed as 


-5 _+x 
at n-1 


Next, the corrected values at time n are determined from 


Ss = E eS 12122 
x 2 1 2 п n 1? 


where Y is the averaging coefficient. Finally, S. is stored in place 
of S. and the model continues on to the next step. 

The value of Y must be carefully chosen as it does effect the com- 
putational stability. As y increases, the maximum time step decreases. 
Haltiner and Williams prefer a Y less than 0.25 in order to permit a 
reasonable time step. Additionally, it should be noted that a large 
value of Y will begin to dampen waves in the meteorological frequency 
range. On the other hand, Robert used a value of y = 0.01 applied over 
a total forecast length of many days. This very weak filter was all 
that was needed for noise suppression. 

For the model in this investigation Y is an input parameter. For 
One set of tests in Chapter X, the value of y is varied from 0.01 to 


0.04. 
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р. DETERMINATION OF TIME STEP 


Cullen (1973) calculated the stability criterion for the two-dimen- 


sional problem as: 


DID. re VI-3 


the minimum grid distance, Ax, varies from 200 km to about 40 km. The 
propagation speed, c, varies from 40 m/sec to 0. Accordingly, the time 
step has been chosen as 60 minutes for the coarse mesh versions, and 


ten or five minutes in the finer mesh versions. 
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VII. NODAL GEOMETRIES 


Once the decision has been made to use triangular elements, the next 
choice is the type of triangle. This investigation concerns both right 


and equilateral triangles. 


A. RIGHT TRIANGLES 

Figure 10 is an example of the right triangle arrangement used in 
this model for a uniform grid. The domain is divided into a series of 
rectangles, and then each rectangle is bisected by a diagonal to form a 
pair of triangles. Kelly and Williams found that significant bias was 
introduced into the model if all of the diagonals sloped in the same 
direction. This model overcomes that kind of bias by alternating the 
slope of the diagonals from rectangle to rectangle both horizontally and 
vertically. Consequently, the connectivity of the nodes will vary from 
four to nine total connections, and each node connects with itself. 

Two methods are used to vary the grid resolutions with right 
triangles. In the first method a FORTRAN DATA statement specifies co- 
efficients used in the calculation of the nodal coordinates. This method 
is very simple and has the advantage of keeping all the nodes in a 
rectangular arrangement.  Rectangularly arranged nodes lend themselves 
well to harmonic analysis with some interpolation. The grid may be 
varied in either one or both the x and y-directions. Figure ll is an 


example of the nodal arrangement by this method with variation in both 
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Figure 10. The domain divided into right triangular elements. Only 
the elements in the lower-left corner are drawn, but the 


pattern is repetitive over the entire domain. Each node 
is marked with an 'X'. 





Figure ll. Same as figure 10 with non-uniform elements. 
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directions. The disadvantage here is that the resolution must be kept 
relatively fine in parts of the domain away from the area of interest 
(the center). Optimumly, the resolution should be fine about the area 
of interest and more coarse elsewhere, hence the next method. 

In the second method, the domain is divided into quadrants with the 
Origin at the center of the domain. The nodes around the outer boundary 
of the first quadrant are specified as having the same coordinates as 
the uniform grid. But the nodes along the abscissa and ordinate are 
compressed toward the origin by the use of another DATA statement.  Con- 
sider four nodes each with coordinates (x. Yi); і = 1...4 . М№ае 1 is 
on the ordinate, node 2 on the right boundary, node 3 on the abscissa, 


and node 4 on the top boundary. The line from node 1 to 2 is defined 


by 


a; 22322122 
i 2 JE 


And the line connecting nodes 3 and node 4 is 


у - у y, -y 
3 
23 4 3 


Equations VII-1 and VII-2 are then solved simultaneously to calculate the 
coordinates of the node located at the intersection of the two lines. By 
selecting the proper pairs of nodes along the boundaries and axes, the 


coordinates of all the nodes on the interior of the first quadrant are 
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found by VII-1 and VII-2. The nodes about the origin are then slightly 
adjusted to allow local uniformity. Finally the first quadrant is re- 
flected about the axes to form the other three quadrants in the domain. 
Figure 12 is an example of such an arrangement. 

The advantage of this method is that the resolution is fine only in 
the area of interest, and changes smoothly in all directions away from 
that area. The disadvantages are the complexity and the non-rectangu- 


larity of the nodes. 


B. EQUILATERAL TRIANGLES 

Preliminary results of the model indicated the presence of a degree 
of noise when uSing right triangles. Cullen (personal conversation) 
suggested the use of equilateral triangles as an alternative element 
arrangement. Figure 13 is an example of such an arrangement. The 
right triangles along the sides of the domain are in reality halves of 
equilateral triangles. Because of cyclic continuity, the domain actually 
"wraps around" so that all of the elements are equilateral. Note that 
the maximum connectivity is now only seven and that the domain is no 
longer square but rectangular. 

Variation of resolution is now somewhat more complicated and is 
accomplished by a transformation of coordinates. The desired end is a 


smooth and easily controllable stretchage and shrinkage of the coordi- 


nate axes. In the x-direction, the transformation is 


Х = х + а соѕ kx VII-3 
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Figure 12. An alternate method of varying resolution. 
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Figure 13. Same as figure 10 except with equilateral triangular 
elements. 


44 \) 





where X = transformed coordinate 
x = original coordinate 


a = coefficient to be determined 


гт 
L 


L = channel length 


X. Р 
the map factor is defined by 


9x = ] - ak sin kx 


Ox 


whose maximum and minimum values are 


IX _ ‚ ‚ӘХ _ _ 
С. poa x min | Ёс 


The ratio, 21, of maximum stretch to minimum shrink is 


= l + ak 
l 1 - ak 


ШИ solving for a 


a= К(кү +1) VII-4 


TO keep the transformed coordinate system from folding back on itself, 


a must be positive. That is, rj, must be equal to or greater than one. 
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Variation of resolution is accomplished by selecting an appropriate 
value for ri and then substituting equation VII-4 into VII-3. Similar 
expressions can be derived for the transformation of the y-coordinate 
by the use of the ratio PE 

The advantage of this method is the extremely sensitive control upon 
resolution afforded by the ratios гу апа r5. Once programmed, this method 
is much easier to use than the more cumbersome DATA statements of the 
previous methods. The major disadvantage is that the transformed system 


is no longer made up of equilateral triangles if one of the ratios is 


different than one. Figure 14 is the nodal arrangement with ку = 5 and 
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Figure 15. Initial S field for the diffusion tests. 
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Vitis | ORCING TERM 


As part of an effort to improve the testing of a numerical scheme, 
it is often beneficial to force the scheme with the desired solution, 
if known. Since this model is a simple advection scheme, the output will 
look very similar to the initial condition. 


Assume the solution is of the form 


S = сов (М(Х = At) + ô] 
where 
2T 

ROS Y 
A = phase speed, 

ТЇ 
ст 
Х = transformation of x coordinate. 


ЇЇ у = 0, then the forced advection equation is 


95 oS _ 

scan ux Ext) VIII-l 
but 

oS _ TT 

DUE uA siníu(X - At) - = VIII-2 
and : 

3S _, ӘХ, T К 

us -u Эх Sin (U(X - At) 5] VIII-3 
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ENSeTEHuting VIII-2 and VIII-3 into VIII-1 


H(A - u 2 sin[u(X - At) - 


3m К = Е(х,С) 


IE 


which is equivalent to 
u(A - u 235) [sin (ux - 7] cos (uàt) - cos (uX - 2) sin(Aut)] = F(x,t) 


applying Galerkin's technique this equation eventually reduces to: 
= g.) <V. ‚> = F vIII-4 
= вар 3! Vi (x€) 


where 


u OX | Т 
E = u(A u Ax) sin (ux 5) cos рле 
VIII-5 
ðX T 
2 = H(A - u Эх” cos (UX =>} Sume 


Now, VIII-4 should be included on the r.h.s. of the predictor equation 
III-7. Notice that both of the equations VIII-5 contain a time depen- 
dent part. The technique now is to calculate the time independent 
coefficients of VIII-5 at the onset and store them for the duration of 
program execution. During the forecast sequence, at each time step, 
the time dependent parts of VIII-5 are calculated, multiplied by the 
applicable pre-stored coefficients and then assembled into the mass 


Matrix as indicated by VIII-4. 
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This forcing term will not be used in most of the model runs as a 
standard part of the predictor equation, but will be the subject of a 
specific set of tests designed to determine its contribution to 


forecast accuracy. 
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IX. ANALYTIC SOLUTION 


Consider the one dimensional equivalent of equation III-l with no 


diffusion: 


3s 3s _ 
3c + u(x) == = 0 IX-1 


where flow is zonal only and the velocity may vary spatially. Recall 


from Chapter V that the zonal wind equation is 


— ТІ 
u(x) = U + B cos(kx - 5) IX-2 
ӘТ 
where k = E '/ L = zonal wave length. 


In a simple one dimensional advection equation for any time t, the 


following holds: 
S(x,t) = S (xav O) 


where х, is some point upstream in the initial field. Now if x, can be 
determined as a function of the current location and the elapsed time, 
then the analytic solution for the forecast at that time would be 


achieved. That is, 


x. = F(x,t) 


S 





Since u = then from IX-2 


ax 
dt ' 


——— sos + ана. — 
d = U B cos (kx, 2) 


Integrating 


kdx 
2 = kdt 


ПП +В (К 2 
COS 9 2 


With the help of a Table of Integrals, this expression is evaluated as 


-— L T 
(U-B)tan 5 (kx pa) 


2 Sl 
{tan I 3 
mar. Tare 
2р (9-8) кап 5(kx, - I) 
Mua 6 
(07 - B`) 
Let 
es U-B 
А 
So that 
1 1 7... -l Шэнэ 1 2-21/2 
сап [x tan >(kx 2)! сап [r tan (kx, 2)1 = „u B ) E 
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Solving for x 


T 1 -- 172 
un 


- ZI J 
x. = — + Z tan (= tan [tan = (= Сап 5(кх- 5) = 2 (9 kt] } IX-3 


2 
The analytic solution at time t for any node at coordinates (x,y) 


0: Тһеп х, апа 


is found by substituting x and t into ТХ-3 to obtain x 
y are substituted into equation V-l. This process requires each value 
of x to be operated upon by five trigonometric functions, each of which 
further compounds truncation error. As a result, this subroutine must 
be run in double precision. Even so, there is still some small error 
in the analytic solution field. This investigation is concerned with 
the relative errors of various nodal arrangements and is not intended 
for comparison studies between finite element and other methods. Since 


the errors in all test cases are relative to the same solution, the 


small error in the solution is not considered significant. 
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Х, RESULTS 


The model contains eight variables: element type (equilateral and 
right triangles), east-west resolution, north-south resolution, wave 
number, diffusion, Robert filter, uniformity of flow, and number of nodes 
on the domain. Additionally, the time step may be varied and the 
initial conditions may be changed so that sharper waves are created. 

If all of these variables were tested throughout their meteorological 
range, and if all the interactions between the variables were investi- 
gated, the number of computer runs required would be well over ten 
thousand. Since that many runs is not practical, decisions have to be 
made concerning the scope of this investigation. There seems to be one 
of two ways to proceed: 

1. One or two questions could be answered rather exactly with an 
attempt to quantify the relationship between only a couple of variables. 
2. General answers can be attempted for many more questions with 
the aim to set qualitative guidelines for future finite element studies. 

In the first alternative above, quantitative relationships would be 
of little value unless they could be extended to FEM models as a class. 
Proof that such an extension is valid is clearly beyond this investiga- 
tion. As a result, the second alternative is chosen and only rather 
general questions are qualitatively answered. 

The model is run a total of 137 times. The runs are compared and 


contrasted in a number of test cases. Several test cases comprise the 
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investigation of each question. Two types of statistics are calculated 
for each run. First of all, a harmonic analysis is accomplished using 
nodes along each latitude circle, for the initial condition, the forecast 
field and the analytic solution field. Since the analysis requires 
equally spaced values, there is an interpolator that calculates values 

at specific points. The linearity of the basis functions lends itself 
very well to linear interpolation. One subroutine determines in which 
element any particular interpolation point lies. Another routine uses 
the coordinates of that point and the values at the nodes of the proper 
element to calculate the interpolated value. In the runs that include 
the forcing term (see Chapter VII), the interpolater is disconnected and 
the raw nodal values are used. This harmonic analysis generated ampli- 
tude and phase information for all possible wave numbers. In the follow- 
ing paragraphs the term"phase speed" denotes the ratio of the phase 

shift of the forecast wave to the phase shift of the wave in the analytic 
solution. 

The second type of statistics generated considers the forecast and 
solution values as an ordered pair for each node. А root-mean-square- 
error (RMSE), correlation and bias are calculated for this set of pairs. 
By-products of these statistics that are not referred to hereafter are 
the applicable means and standard deviations. 

The following sections are organized such that each one addresses 
an individual problem, question or technique. To this end, within each 
Section, one or more of the eight variables listed above is varied from 


a standard or nominal model configuration. Changes in the generated 
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meeaersctics are then attributed to the influence of the altered variables. 
The nominal model configuration follows: equilateral elements, uniform 


= 1, see Chapter VII), wave number one, no diffusion, 


grid (ri =r, 


Robert filter, y = 0.1, wind perturbation ratio R = 0.0 (see Chapter V), 
number of nodes = 13 Х 12. Тһе standard forecast length is 48 hours, 
during that period any feature, moving with the mean flow, should move 
3456 km or around the domain 1.44 times. Section A is a general statement 
based upon nearly all of the runs. It includes changes made to all of 

the variables. In all of the remaining sections, the departure from the 


nominal configuration is annotated. 


A. ACCURACY VS. RESOLUTION 

For this test, six runs are discarded because they are deliberately 
constructed to test unique factors discussed later. In the other 131 
runs, the forecast values stay within reasonable limits for all combina- 
tions of resolution. In general, the accuracy over the domain as a whole 
decreases as the stretch of the grid resolution increases, but the 
statistics all remain bounded. The correlations are all above 0.96, the 
phase speeds are within four percent and the RMSE is always less than 

eight percent of the range of the amplitude. Typical values are cor- 
relation 0.99, phase speed 1.008 and RMSE one percent of range. 

This result is very encouraging in that no cases have to be dis- 
regarded simply because they can not be explained. The cases are 


consistent with one another and exceptions are few. 
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БО EFFECT OF DIFFUSION 

A particularly noisy case is used for this test. The domain length 
scale has been doubled so that the non-uniform wind is more obvious. 
With resolution of г, = 1, r, = 2 wind perturbation ratio R = 0.4, and 
13 X 24 nodes, the diffusion values are incrementally increased from 
zero. A total of eight runs are made. Figure 15 is the initial condi- 
tion for this test and Figure 16 is the 48 hour forecast with no diffu- 
sion. As diffusion is increased the reduction in RMSE is dramatic and 
the correlation increases, but the wave amplitude diminishes and the 
maximum time step must be decreased to insure stability. However, a 
critical point is eventually reached at which the amplitude is reduced 
so much that the RMSE starts to climb again and correlation falls. 
Figure 17 is the 48 hour forecast with the value of diffusion that 
nearly minimizes the RMSE. Note the absence of noise but the obvious 
loss of amplitude. Compare this with Figure 18, the analytic solution. 
Notice that the contour packing is downstream of the high in Figures 16 
and 18, but upstream of the high in Figure 17. This reversal of packing 
is a by-product of the manner in which the model applies diffusion. 
Cullen (personal communication) suggests that diffusion be increased in 
the areas of fastest flow and decreased where the flow is slower. 
Following this advice, the model applies diffusion proportional to the 
windspeed at each node. The effect with non-uniform flow is that diffu- 
Sion is also non-uniform, hence the slight alteration in the forecast 
field. Diffusion appears to be a very effective noise filter, but tests 
must be made independently for each configuration because the useful 


limit of diffusion will vary with nodal geometry. 
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ENwure l6. 48 hour forecast for field in figure 15 with no diffusion. 
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bagure 17. 


Same as figure 16 but with optimum diffusion. 
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Figure 18. Analytic solution for diffusion tests. 
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Figure 19. Initial S field for Robert filter tests. 
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fee EFFECT OF ROBERT FILTER 

For this test the model is configured nominally to maximize the accu- 
racy of the solution. The only variable that is changed is the Robert 
averaging coefficient. Figure 19 is this initial configuration. The 
averaging coefficient, Y starts at a value of 0.4 and is slowly decreased 
to 0.01 over ll runs. The RMSE decreases from 8.7 to 1.74 over this range 
with a corresponding increase in correlation and the phase speed approaches 
unity. 

Since the filter is desirable to control high frequency noise, a com- 
promise must be reached that will not unduly degrade the forecast. For 
the rest of this investigation a coefficient of y = 0.1 is used since the 
associated RMSE is only 2.54. Figure 20 is the 48 hour forecast field 
with this coefficient and Figure 21 is the corresponding analytic solution. 
There is little improvement in RMSE left to be realized if any effective 


filtering is desired. 


D. ACCURACY VS. WAVE NUMBER 

Again the model is configured nominally, except that the number of 
nodes is 13 X 24. Six runs are made with the wave number varying from 
one to eight. One would expect the forecast to deteriorate as the wave 
number increases. This is due in part to the fewer number of nodes along 
a latitude circle available to resolve each wave. But the deterioration 
of RMSE is primarily the result of the displacement error measured in 
terms of the wavelength. A small displacement error with a short wave 
would produce a much greater deviation from the true solution than would 


be produced by the same displacement error in a long wave. 
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48 hour forecast for field in figure 19. 


Figure 20. 


Analytic solution for Robert filter tests. 





Figure 21. 
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The results of this test verifies this reasoning. As the wave number 
of the initial field increases the RMSE worsens, even though the phase 
speed error stays within four percent. Four of these cases were excluded 


from Section A above. 


Phase Speed 
Wave Number Ratio 


1.008 
12008 
1202 
17203 
1.04 
103 


Om & WN KF 


E. RIGHT TRIANGLES VS. EQUILATERAL TRIANGLES 

Five pairs of runs are used to test the hypothesis that elements 
formed by equilateral triangles yield forecasts superior to those based 
upon elements formed by right triangles. In all cases the flow is 
uniform but the resolution is varied among the cases from uniformity 
to variations in either or both directions. Partial results appear in 
Table I as Cases 1 through 5. Figures 22, 23 and 24 are the initial 
field, 48 hour forecast and analytic solution, respectively for the 
right triangles in Case 3. Figures 25, 26 and 27 are the corresponding 
figures for equilateral triangles also in Case 3. 

The differences between the two arrangements, while not dramatic, 
are significant and consistent. In every case the RMSE is lower and 
correlation higher for equilateral triangles. The average reduction in 


RMSE is about 20 percent. Not shown in Table I but readily apparent in 


62 





2400 


km 





0 2400 


Figure 22. Initial S field for Case 3 with right triangles. 
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Figure 23. 48 hour forecast for field in figure 22. 
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Analytic solution corresponding to figure 23 


Figure 24. 


Initial S field for Case 3 with equilateral triangles. 





Figure 25. 
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Figure 26. 


48 hour forecast for field in figure 25 
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Figure 27. 


Analytic solution corresponding to figure 26. 
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the harmonic analysis, is the generation of a spurious wave number five 


in the right triangle arrangement in four of the five cases. 


TABLE I 


RESULTZSTOFZRIGHT VS. EZEOUTEATERAL ELEMENT 
AND DIRECTION OF FLOW TESTS 


RIGHT EQUILATERAL 
E E 

CASE 1 2 RMSE CORR. RMSE CORR. 
ШІ 1 1 4,25 .9994 2.54 ,9999 
2 2 iL 4.54 .9993 4.17 .9996 
3 1 2 6.40 .9991 4.87 .9996 
4 2 2 6.96 .9990 5.60 .9994 
5 3.55 1 6.19 .9981 5.75 .9992 
6 i 4 7.48 .9980 
7 4 4 6.60 .9990 


8 1! 1552 279923950 


F. EFFECT OF ADDING NODES 

Next, 16 runs are combined to form nine test cases in order to 
determine the effect that adding nodes has on the accuracy of the fore- 
cast. In this section, the variables are wave number, resolution and 
number of nodes. It is found that where the model does very well (long 
waves, uniform flow, relatively small difference in resolution over the 
entire domain), the addition of more nodes does not increase the accuracy, 
especially if the additional nodes destroy the original equilateral 
triangularity. However, if the nodes are added symmetrically in both 
the north-south and the east-west direction, so that the effect is to 
increase the resolution uniformly while preserving equilateral triangu- 


larity, then the result is an improved forecast. 
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On the other hand, when the mndel is initialized with shorter waves, 
more nodes will improve the forecast even though they may not be added 
symmetrically. The major advantage of additional nodes is that they 
allow the resolution to be changed more slowly and more smoothly while 
attaining a desired minimum resolution. The arrangement with fewer nodes 
requires a faster and more abrupt resolution change in order to attain 
the same minimum resolution. The smoothness of the transition between 
maximum and minimum resolution is apparently the critical factor. Later 
Paragraphs also address this factor. The phase speed ratios for the 


uniform grids are as follows: 


Wave Phase Speed 
Nodes Number Ratio RMSE 
23412 1 1.001 21: 
132x924 1 1.008 10.4 
259% 24: 1 1.003 7.4 
19322212 2 15003 60.29 
13 x 24 2 12002 30.7 
25 356 24 2 1.003 2272 


G. VARIATION OF RESOLUTION VS. DIRECTION OF FLOW 

The initial conditions for the model are deliberately set to allow 
only zonal flow. (See Chapter V.) One of the reasons for this is to 
be able to specifically test the effects of variation of resolution along 
the direction of flow. А total of 13 runs form six test cases. Many of 
these runs are the same ones as used in tests in preceding paragraphs. 
Some of these cases appear in Table I. In that table the comparisons 
are now vertical. For example, compare Cases 2 with Case 3 to determine 


that ri = 2 and г. = 11S superior to r. = 1 апа r, = 2.106 bobh right 
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and equilateral triangles. This means that varying the resolution with 
the flow is superior to varying across it. Case 2 is illustrated in 
Figures 28 and 29 which are the initial field and 48 hour forecast re- 
spectively. They should be compared with Figures 25 and 26 which, as 
stated previously, are part of Case 3. The same result holds for every 
case on the table. Additional tests have also been run with resolution 
ratios of six and eight. The results are the same. 


One further test can be conducted with the data on Table I. Resolu- 


оп ratios of r, = r 2 will yield the same overall resolution as one 


1 2 
of the ratios equal to 3.55 and the other unity. This is the rational 
for Cases 4, 5 and 8. These cases show that better results can be ob- 
tained by varying resolution in both directions as opposed to only across 
the flow by an equivalent amount. This fact agrees with the results of 
ehe previous section in that varying resolution in both directions tends 


to maintain the equilateral triangularity better than varying across the 


flow alone. 


H. SMOOTH VS. ABRUPT VARIATION OF RESOLUTION 

For this test, a version of the model is employed that allows the 
resolution to be changed abruptly. That version is paired with a more 
Standard version that allows smoother variation of resolution. Both 
versions use right triangular elements. Both versions generate the same 
Maximum resolution, but with very different transitions to the area of 
coarser mesh. Twelve runs are paired into six tests. Two different 


Wave numbers and three kinds of flow are tested. These runs included 
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Figure 28. 
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48 hour forecast for Case 2 


Figure 29. 
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both the normal and "flipped" (See Section J below) configuration. The 
results are contained in Table II. Even though only six tests are used, 
the results are dramatic and highly significant. The reduction in RMSE 
from the abrupt to the smooth cases averages 62 percent. Although not 
shown, the reduction in bias averages 70 percent. Additionally the 
abrupt cases generate significantly more noise at all frequencies. Cases 
10 and 12 are the last two cases excluded from Section A above. Figures 
30 through 33 illustrate the evolution of the forecast in Case 9. The 
figures are the initial field and 12, 24 and 48 hour forecast respective- 
ly for the abrupt change. Compare these to Figures 34 and 35 which are 
the initial field and 48 hour forecast for the corresponding smooth change. 

Combining this result with that of Section F above, it is obvious 
that the model is highly sensitive to the rate and smoothness of the 


eransıtion from fine to coarse resolution. 


TABLE II 


RESULTS OF SMOOTH VS. ABRUPT TESTS 


WAVE WIND SMOOTH ABRUPT 

CASE NUMBER RATIO RMSE CORRE RMSE CORR. 
9 1 0.0 4.54 7999 13,5 1989 
lO 2 0.0 283 .956 S352 ‚470 
iE 1 2 9.46 2925 29,5 . 968 
12 1 0.4 21.4 2979 51572, „895 
ШЕ 1 022 4.53 299 14.6 227511) 
ШОР 1 0.4 8-83 2995 20.0 2979 


Т. EFFECT OF VARIABLE WINDS 
In several of the previous tests, reference has been made to non- 


uniform flow. For this section and the following two sections a total 
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Figure 30. Initial S field with abrupt change of resolution 
(Case 9). 





Figure 31. 12 hour forecast for field in figure 30. 
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EUgure 32. 


48 hour forecast for field in figure 30. 
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Figure 33. 








Figure 34. Initial S field with smooth change of resolution 
(Case 9). 
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Figure 35. 48 hour forecast for field in figure 34. 
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ОЁ 42 runs are made to study the effects of wind and the forcing term. 
The 42 runs are formed into 14 tests in which the wind ratios are 
specified as 0.0 for uniform flow and 0.2 and 0.4 for variable flow. 
The tests include cases that have both uniform and non-uniform resolu- 
tion, both equilateral and right triangular elements and the 13 x 12 
and 13 x 24 nodal arrangements. Table III contains the results of 12 
of the runs which comprise four tests. In 12 of the total of 14 tests, 
the model performs better with the uniform wind. In all tests once 
the wind is non-uniform the forecasts are degraded as non-uniformity 


increases, 


TABLE III 


RESULTS OF VARIABLE WIND, "FLIPPED" WIND AND FORCING TERM TESTS 


UNEORCED FORCED 
CASE RATIO RMSE SPEED CORR. RMSE SPEED CORR. 
NEED... NORMAL oa aaa ©.0l e166 G6 ee 6 ee ke eee 0.0 6 eed wee o e686 
ES 0.0 21:28, 1.004 49996 247 1.004 225996 
14 022 3.4 1559005 19905 2759 15095 59945 
15 0.4 957 122121215: 22970 eM 1.006 29970 
Ниэттэгэээ ОИ. БЕТЕРЕР о СОНИ о оо оо сое ое ео о оо ооо ооо оо во 
16 0.0 268 1.003 299099 2.8 1:007 23996 
17 12 Ze 1.004 „9996 276 1.004 (9997 
18 0,4 S 1. 005 9935 2-3 1,004 29960 


J. "FLIPPED" WINDS 

Up until now, the initial conditions, in conjunction with the trans- 
formation of coordinates scheme, require that the area of highest wind 
is also the area of highest resolution, provided, of course, that both 
the wind and the resolution are non-uniform. Now, downstream of the area 


of highest winds 1s an area of convergence and downstream of the area 
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of lowest wind is an area of divergence. Therefore, where the conver- 
gence runs into the divergence is the area where the gradients are 
highest. It is a "bottleneck". But that is exactly the area where the 
resolution is the most coarse. The forecast should improve if the grid 
or initial wind pattern could be "flipped" to concentrate the nodes in 
the area of minimum winds. Since "flipping" the grid requires fewer 
changes, that is the choice here. As an example, Figures 36 and 37 are 
the initial field and 48 hour forecast of such an arrangement. 

To verify this reasoning, 24 runs are formed into 12 tests. Cases 
13, 14, and 15 on Table III are "normal" whereas Cases 16, 17 and 18 
are "flipped". In all tests, the "flipped" arrangement yields superior 
forecasts. This result justifies the reasoning that the high resolution 
should be concentrated in area of highest gradients, not the areas of 


fastest flow. 


КОО FORCING TERM 

The forcing term was developed in Chapter VIII as an alternate method 
of testing the forecast if the general form of the end product is known. 
The 42 runs form 21 pairs of forced/unforced tests, some of which are in 
Table III. The forced forecasts are slightly better in 15 tests, the 
unforced are slightly better in five tests. One test comes out a tie. 
In general, the unforced forecasts are better if the wind is uniform while 
a non-uniform wind favors the forced forecast. However, all differences 
are very slight. Figures 38 and 39 are the initial field and 48 hour fore- 
cast respectively for a typical forced forecast with a uniform wind. They 
should be compared with Figures 28 and 29 which form the corresponding 


unforced case, 


75 





© 
^ 
o 
N 





E 
5 
ж 
о 
Figure 36. Initial S field with "flipped" winds. r. = 2595 
ib 2 
and R = 0.2. 


= 336 
32:15! 
y WU E 





Ergure 37. 


48 hour forecast for field in figure 36. 
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Initial S field for a typical forced case. 


Ergure 38. 





48 hour forecast for field in figure 38 


Figure 39. 
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XI. CONCLUSIONS 


This report applies the Finite Element Method to the two-dimensional 
advection equation with diffusion. Although the underlying equation is 
Simple enough, FEM models have certain inherent advantages, which include 
very accurate phase speeds and the ability to change resolution easily. 
However, they also have limitations, such as a relatively high computa- 
tional cost, and a requirement for fairly sophisticated and efficient 
programming. The purpose of this investigation is to establish some guide- 
lines and provide qualitative methods that may be used in the development 
of future FEM models. 

The conclusions are given below in the same order as they are discussed 
in Chapter X: 

l. The model yields reasonable and consistent results. There appear to 
be no major algorithmic or coding problems as all the results are 
mathematically and meteorologically plausible. 

2. Diffusion is an effective noise filter; however it does have a maximum 
usable limit which must be independently determined for each applica- 
tion. Numerical models should contain at least some provision for 
diffusion. 

3. The Robert filter is also a valuable feature to have built into a 
model. In general, the filter should be no heavier than absolutely 
needed to control high frequency noise. This model uses a filter of 


0.1; higher values are not recommended. 
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10. 


ШТ, 


The model handles wave numbers l and 2 well (24 and 12 nodes per 
wave, respectively) and wave number 3 (8 nodes per wave) with less 
accuracy. Higher wave numbers maintain accurate phase speeds 
although RMSE deteriorates. 

Formulations with elements as equilateral triangles outperform those 
with elements as right triangles. The difference in RMSE is about 
20 percent. The right triangle formulation generates spurious short 
waves which may be be the primary cause of the relative inaccuracy. 
Adding nodes (and hence, elements) to the domain will improve the 
model if they are added symmetrically. More nodes allow smoother 
transition between fine and coarse resolution. This smoothness is 
erleical. 

An unexpected result is that this model performs better when the reso- 
lution is varied with the flow rather than across it. The reasons 
for this are unclear. 

The most important conclusion of this investigation is that the 
smoother and slower the change in the resolution, the better the 
forecast. In this model the smooth change reduces the RMSE by 62 
percent over the abrupt change. 

The model works best with uniform flow. In general, as the flow 
departs from uniformity, the forecast continues to degrade. 

The forecast can be improved by concentrating the high resolution 
areas in the region of the strongest gradients. 

Inclusion of a term that forces the supposed solution back into the 
model may yield slightly better forecasts. It forms a valid, alterna- 


tive method to test numerical models. 
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